library(tidyverse)
library(janitor)
library(xtable)
library(broom)
library(lfe)
library(stargazer)

lics <- read_rds("glick_palmer_ozs.rds")

tabyl(lics, lic_type)
tabyl(lics, lic_type, qoz)
tabyl(lics, lic_type, qoz) %>% adorn_percentages()

lics <- filter(lics, lic_type!="contig")
lics <- mutate(lics, gov_party=ifelse(gov_rep==1, "Republican", "Democrat"))
lics <- filter(lics, pop>0)

# Table A1 ----------------------------------------------------------------
x <- lics %>% filter(st!="AK", st!="WV") %>%
  select(gov_vote_cty,ld_party_match_any,lics_in_county_inv,med_income_hh,poverty_rate,ui_investment_score) %>%
  summarize(across(everything(), list(mean= ~ mean(., na.rm=T),
                                      sd=~sd(., na.rm=T),
                                      min=~min(., na.rm=T),
                                      max=~ max(., na.rm=T),
                                      n=~sum(!is.na(.))),
                   .names="{.col}.{.fn}")) %>%
  pivot_longer(everything(), names_sep="\\.", names_to=c("var", "stat")) %>%
  pivot_wider(names_from=stat, values_from=value)
x$var <- c("Gov County Vote %",
           "LD Party-Match",
           "LICs In County (Inv)",
           "Median HH Income",
           "Poverty Rate",
           "UI Investment Score")  
x %>% xtable() %>%
  print(include.rownames=F, only.contents=T, 
        file="tables/table_a1.tex")


# Figure A1 ---------------------------------------------------------------
lics %>% filter(st!="AK", st!="WV") %>%
  group_by(st) %>%
  summarize(over50=mean(gov_vote_cty>.50),
            n=n()) %>%
  arrange(over50, desc(st)) %>%
  mutate(st=factor(st, levels=st)) %>%
  ggplot(aes(y=st, x=over50)) +
  geom_vline(xintercept=.25, color="red") +
  geom_point() +
  labs(y="", x="Percent of LICs in Counties Governor Won") +
  scale_x_continuous(limits=c(0,1), labels=scales::percent_format())
ggsave("figures/fig_a1.pdf", width=5, height=6, units="in")


# Figure 1 / Table A2 -----------------------------------------------------
x <- lics %>% filter(st!="AK", st!="WV") %>%
  select(st, qoz, gov_vote_cty, ld_party_match_any,
         lics_in_county_inv, ui_investment_score, 
         med_income_hh, poverty_rate) %>%
  pivot_longer(c(-st, -qoz), names_to="term", values_drop_na = T) %>%
  group_by(st, term) %>%
  filter(! (st=="RI" & term=="ld_party_match_any")) %>%
  do(tidy(t.test(value ~ qoz, data=.)))

df2 <- x %>% mutate(r=ifelse(p.value<.05 & estimate>0, -2.5, 
                             ifelse(p.value<.05 & estimate<0, 2.5, 0))) %>%
  arrange(term, r, st) %>%
  group_by(term, r) %>%
  mutate(n=n(),
         gr=ifelse(n %% 2 == 0, 
                   ifelse(row_number() %% 2==1, -1, 1),
                   ifelse(row_number()==1, -1, ifelse(row_number() %% 2==1, 1, -1))),
         gr=ifelse(n==1 & r==-2.5, 1, gr)) %>%
  group_by(term, r, gr) %>%
  mutate(y=n()-row_number()+.5,
         x=r+.5*gr)

gov_parties <- select(lics, st, gov_party) %>% distinct()
df2 <- left_join(df2, gov_parties)
df2 <- df2 %>% ungroup() %>%
  mutate(term=factor(term,
                     levels=c("gov_vote_cty",
                              "ld_party_match_any",
                              "lics_in_county_inv",
                              "med_income_hh",
                              "poverty_rate",
                              "ui_investment_score"),
                     labels=c("Gov County Vote %",
                              "LD Party-Match",
                              "LICs In County (Inv)",
                              "Median HH Income",
                              "Poverty Rate",
                              "UI Investment Score")))
ggplot(df2, aes(x=x, y=y, label=st, fill=gov_party)) +
  geom_tile(width=1, height=1, color="white") +
  geom_text(color="white", size=2) +
  scale_x_continuous(limits=c(-3.5, 3.5), breaks=c(-2.5, 0, 2.5), 
                     labels=c("Negative", 
                              "None", 
                              "Positive")) +
  scale_fill_manual(values=c("gray10", "gray50"),
                    name="", labels=c("Democratic Governor", "Republican Governor")) +   #c(pal3[1], pal3[3])) +
  facet_wrap(~term, scales="free_x") +
  theme_minimal() +
  theme(panel.grid=element_blank(),
        panel.border=element_rect(color="black", fill=NA),
        axis.text.y=element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x=element_text(size=6, color="black"),
        legend.position = "bottom",
        legend.key.size = unit(.75, "lines"),
        legend.text = element_text(size=6))
ggsave("figures/fig_1.pdf", device=cairo_pdf, width=6,
       height=6, units="in")

x %>% mutate(s = format(round(estimate,2), digits=2),
             s = ifelse(p.value<.05, paste0(s, "*"), s),
             s = ifelse(p.value<.01, paste0(s, "*"), s)) %>%
  select(st, term, s) %>%
  pivot_wider(names_from=term, values_from=s) %>%
  xtable() %>%
  print(include.rownames=F, only.contents=T, 
        file="tables/table_a2.tex")




# Table 1 -----------------------------------------------------------------

m1 <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
           data=lics)
m2 <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
           data=filter(lics, gov_rep==1))
m3 <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
           data=filter(lics, gov_rep==0))

stargazer(m1, m2, m3,
          keep.stat = c("n","rsq","adj.rsq"),
          column.labels = c("All", "Republicans", "Democrats"),
          dep.var.caption = "",
          dep.var.labels = "",
          dep.var.labels.include = F,
          omit="Constant",
          covariate.labels = c(
            "Gov. Vote County",
            "LD Party Match",
            "LICS in County Inv.",
            "Med HH Income",
            "Poverty Rate",
            "UI Investment Score"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001", "Models include state FEs.","Standard errors clustered by state."), 
          notes.append = F, 
          out = "tables/table_1.tex",
          float = F
)


# Table A9 ----------------------------------------------------------------
m1 <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate | state | 0 | state, 
           data=lics)
m2 <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate | state | 0 | state, 
           data=filter(lics, gov_rep==1))
m3 <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate | state | 0 | state, 
           data=filter(lics, gov_rep==0))

stargazer(m1, m2, m3,
          keep.stat = c("n","rsq","adj.rsq"),
          column.labels = c("All", "Republicans", "Democrats"),
          dep.var.caption = "",
          dep.var.labels = "",
          dep.var.labels.include = F,
          omit="Constant",
          covariate.labels = c(
            "Gov. Vote County",
            "LD Party Match",
            "LICS in County Inv.",
            "Med HH Income",
            "Poverty Rate"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001", "Models include state FEs.","Standard errors clustered by state."), 
          notes.append = F, 
          out = "tables/table_a9.tex",
          float = F
)


# Tables A3, A4, A5 -------------------------------------------------------
lics <- mutate(lics, lics_in_county_ln=log(lics_in_county),
               lics_in_county_4 = ifelse(lics_in_county<4, 0, 1))

m1a <- felm(qoz ~ gov_vote_cty + ld_party_match_any + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=lics)
m1b <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=lics)
m1c <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_ln + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=lics)
m1d <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_4 + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=lics)
m1e <- felm(qoz ~ gov_vote_cty + ld_party_match_any + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, st!="TX", st!="CA"))
m1f <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, st!="TX", st!="CA"))
m1g <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_ln + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, st!="TX", st!="CA"))
m1h <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_4 + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, st!="TX", st!="CA"))

m2a <- felm(qoz ~ gov_vote_cty + ld_party_match_any + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==1))
m2b <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==1))
m2c <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_ln + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==1))
m2d <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_4 + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==1))
m2e <- felm(qoz ~ gov_vote_cty + ld_party_match_any + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==1, st!="TX", st!="CA"))
m2f <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==1, st!="TX", st!="CA"))
m2g <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_ln + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==1, st!="TX", st!="CA"))
m2h <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_4 + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==1, st!="TX", st!="CA"))

m3a <- felm(qoz ~ gov_vote_cty + ld_party_match_any + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==0))
m3b <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==0))
m3c <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_ln + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==0))
m3d <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_4 + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==0))
m3e <- felm(qoz ~ gov_vote_cty + ld_party_match_any + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==0, st!="TX", st!="CA"))
m3f <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==0, st!="TX", st!="CA"))
m3g <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_ln + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==0, st!="TX", st!="CA"))
m3h <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_4 + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==0, st!="TX", st!="CA"))

stargazer(m1a, m1b, m1c, m1d, m1e, m1f, m1g, m1h,
          keep.stat = c("n","rsq","adj.rsq"),
          dep.var.caption = "",
          dep.var.labels = "",
          dep.var.labels.include = F,
          omit="Constant",
          covariate.labels = c(
            "Gov. Vote County",
            "LD Party Match",
            "LICs in County",
            "LICs in County ln",
            "LICS in County > 3",
            "Med HH Income",
            "Poverty Rate",
            "UI Investment Score"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001", "Models include state FEs and standard errors clustered by state. Models 5-8 exclude California and Texas"), 
          notes.append = F, 
          out = "tables/table_a3.tex",
          float = F
)

stargazer(m2a, m2b, m2c, m2d, m2e, m2f, m2g, m2h,
          keep.stat = c("n","rsq","adj.rsq"),
          dep.var.caption = "",
          dep.var.labels = "",
          dep.var.labels.include = F,
          omit="Constant",
          covariate.labels = c(
            "Gov. Vote County",
            "LD Party Match",
            "LICs in County",
            "LICs in County ln",
            "LICS in County > 3",
            "Med HH Income",
            "Poverty Rate",
            "UI Investment Score"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001", "Models include state FEs and standard errors clustered by state. Models 5-8 exclude Texas."), 
          notes.append = F, 
          out = "tables/table_a4.tex",
          float = F
)

stargazer(m3a, m3b, m3c, m3d, m3e, m3f, m3g, m3h,
          keep.stat = c("n","rsq","adj.rsq"),
          dep.var.caption = "",
          dep.var.labels = "",
          dep.var.labels.include = F,
          omit="Constant",
          covariate.labels = c(
            "Gov. Vote County",
            "LD Party Match",
            "LICs in County",
            "LICs in County ln",
            "LICS in County > 3",
            "Med HH Income",
            "Poverty Rate",
            "UI Investment Score"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001", "Models include state FEs.","Standard errors clustered by state. Models 5-8 exclude California"), 
          notes.append = F, 
          out = "tables/table_a5.tex",
          float = F
)


# Table A6 ----------------------------------------------------------------

lics %>% filter(st!="AK", st!="WV", st!="NE") %>% 
  group_by(st) %>%
  nest() %>%
  mutate(m=map(data, ~ lm(qoz ~ gov_vote_cty + ld_party_match_any + 
                            lics_in_county_inv + med_income_hh +  
                            poverty_rate + ui_investment_score,
                          data=.) %>% tidy())) %>%
  unnest(cols=m) %>%
  mutate(s = format(round(estimate,3), digits=2),
         s = ifelse(p.value<.05, paste0(s, "*"), s),
         s = ifelse(p.value<.01, paste0(s, "*"), s)) %>%
  select(st, term, s) %>%
  pivot_wider(names_from=term, values_from=s) %>%
  select(-`(Intercept)`) %>%
  xtable() %>%
  print(include.rownames=F, only.contents=T,
        file="tables/table_a6.tex")


# Table A7 ----------------------------------------------------------------

m1s <- felm(qoz ~ gov_vote_swing + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=lics)
m2s <- felm(qoz ~ gov_vote_swing + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==1))
m3s <- felm(qoz ~ gov_vote_swing + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
            data=filter(lics, gov_rep==0))

stargazer(m1s, m2s, m3s,
          keep.stat = c("n","rsq","adj.rsq"),
          column.labels = c("All", "Republicans", "Democrats"),
          dep.var.caption = "",
          dep.var.labels = "",
          dep.var.labels.include = F,
          omit="Constant",
          covariate.labels = c(
            "Gov. Swing County",
            "LD Party Match",
            "LICS in County Inv.",
            "Med HH Income",
            "Poverty Rate",
            "UI Investment Score"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001", "Models include state FEs.","Standard errors clustered by state."), 
          notes.append = F, 
          out = "tables/table_a7.tex",
          float = F
)


# Table A10 ---------------------------------------------------------------

m1 <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate + ui_investment_score + ur_code | state | 0 | state, 
           data=lics)
m2 <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate + ui_investment_score + ur_code | state | 0 | state, 
           data=filter(lics, gov_rep==1))
m3 <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate + ui_investment_score + ur_code | state | 0 | state, 
           data=filter(lics, gov_rep==0))

stargazer(m1, m2, m3,
          keep.stat = c("n","rsq","adj.rsq"),
          column.labels = c("All", "Republicans", "Democrats"),
          dep.var.caption = "",
          dep.var.labels = "",
          dep.var.labels.include = F,
          omit="Constant",
          covariate.labels = c(
            "Gov. Vote County",
            "LD Party Match",
            "LICS in County Inv.",
            "Med HH Income",
            "Poverty Rate",
            "UI Investment Score",
            "UR Class. Large Fringe Metro",
            "UR Class. Medium Metro",
            "UR Class. Small Metro",
            "UR Class. Micropolitan",
            "UR Class. Non-core"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001", "Models include state FEs.","Standard errors clustered by state."), 
          notes.append = F, 
          out = "tables/table_a10.tex",
          float = F
)


# Fig A2 ------------------------------------------------------------------

x <- lics %>% group_by(st, county) %>% summarize(qoz=max(qoz)) %>%
  group_by(st) %>%
  summarize(prop=mean(qoz))
mean(x$prop)
sum(x$prop==1)

x %>% arrange(prop, desc(st)) %>%
  mutate(st=factor(st, levels=st)) %>%
  ggplot(aes(y=st, x=prop)) +
  geom_point() +
  scale_x_continuous(limits=c(0,1), labels=scales::percent_format())+
  labs(y="", x="Proportion of Counties With at Least One QOZ")
ggsave("figures/fig_a2.pdf", width=5, height=6, units="in")


# How many counties had at least one LIC?
actual_num_counties <- lics %>%
  group_by(st, county) %>%
  summarize(qoz=max(qoz)) %>%
  ungroup() %>%
  summarize(qoz=sum(qoz))

x <- lics %>% group_by(st, county) %>% summarize(qoz=max(qoz)) %>%
  group_by(st) %>%
  summarize(prop=mean(qoz))
mean(x$prop)
sum(x$prop==1)



# Table A8 -----------------------------------------------------------------

lics <- read_rds("glick_palmer_ozs.rds")
lics <- mutate(lics, gov_party=ifelse(gov_rep==1, "Republican", "Democrat"))
lics <- filter(lics, pop>0)

m1 <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
           data=lics)
m2 <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
           data=filter(lics, gov_rep==1))
m3 <- felm(qoz ~ gov_vote_cty + ld_party_match_any + lics_in_county_inv + med_income_hh +  poverty_rate + ui_investment_score | state | 0 | state, 
           data=filter(lics, gov_rep==0))

stargazer(m1, m2, m3,
          keep.stat = c("n","rsq","adj.rsq"),
          column.labels = c("All", "Republicans", "Democrats", "Tex.", "N. Car."),
          dep.var.caption = "",
          dep.var.labels = "",
          dep.var.labels.include = F,
          omit="Constant",
          covariate.labels = c(
            "Gov. Vote County",
            "LD Party Match",
            "LICS in County Inv.",
            "Med HH Income",
            "Poverty Rate",
            "UI Investment Score"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001", "Models include state FEs.", "Standard errors clustered by state."), 
          notes.append = F, 
          out = "tables/table_a8.tex",
          float = F
)
